{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kinetic Mechanism Pathway Analyzer\n",
    "\n",
    "This analyzer is designed to provide a better user experience and improved analyzing capacity (e.g., large mechanisms) compared with Chemkin mechanism analyzer. This analyzer takes reaction ROP file (e.g. Chemkin ckcsv file), builds a reaction network, helps answer\n",
    "\n",
    "- if the mechanism has certain pathway (from start species to end species)\n",
    "\n",
    "- what is the biggest flux from/towards a certain species\n",
    "\n",
    "When using this script, please pay attention to the sections marked with **[user input]** which require user to specify certain input parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from rmgpy.rmg.model import CoreEdgeReactionModel\n",
    "from rmgpy.chemkin import loadChemkinFile\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [user input] Prepare files\n",
    "\n",
    "Please add \n",
    "\n",
    "- annotated chemkin files (mechanism and species dictionary) and \n",
    "\n",
    "- chemkin simulation ckcsv file \n",
    "\n",
    "to a folder with name stored in `mech` variable in the next cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mech = 'minimal'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mechPath = os.path.join('./', 'data', 'pathway_analysis', mech)\n",
    "chemkinPath= os.path.join(mechPath, 'chem.inp')\n",
    "dictionaryPath = os.path.join(mechPath, 'species_dictionary.txt')\n",
    "ckcsvPath= mechPath + '/CKSoln.ckcsv'\n",
    "\n",
    "model = CoreEdgeReactionModel()\n",
    "model.core.species, model.core.reactions = loadChemkinFile(chemkinPath,dictionaryPath)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# generate paris for reactions that don't have flux pairs\n",
    "for rxn in model.core.reactions:\n",
    "    if not rxn.pairs: rxn.generatePairs()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create reaction network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import networkx as nx\n",
    "import matplotlib.pyplot as plt\n",
    "from rmgpy.tools.extractInfoFromckcsv import getConcentrationDictFromCKCSV, getROPFromCKCSV, getFluxGraphEdgesDict\n",
    "from rmgpy.chemkin import getSpeciesIdentifier\n",
    "from IPython.display import display\n",
    "import numpy as np\n",
    "from rmgpy.rmg.pdep import PDepReaction\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### step1: prepare a dict for graph edges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "firstColDict, spc_mf_dict = getConcentrationDictFromCKCSV(ckcsvPath)\n",
    "first_col_dict, spc_total_rop_dict, spc_rop_dict = getROPFromCKCSV(ckcsvPath)\n",
    "\n",
    "graph_edges_dict = getFluxGraphEdgesDict(spc_rop_dict, model.core.reactions)\n",
    "\n",
    "graph_edges_dict_simple = {}\n",
    "for pair in graph_edges_dict:\n",
    "    node1 = getSpeciesIdentifier(pair[0])\n",
    "    node2 = getSpeciesIdentifier(pair[1])\n",
    "    graph_edges_dict_simple[(node1, node2)] = graph_edges_dict[pair]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### [user input] step2: pick a time to analyze pathways\n",
    "\n",
    "Please specify \n",
    "\n",
    "- the time point you want to investigate\n",
    "\n",
    "- the chemkin name of your initial species/reactant\n",
    "\n",
    "- the reaction temperature in K\n",
    "\n",
    "- the reaction pressure in Pa if the mechanism includes pressure-dependent reactions\n",
    "\n",
    "- showIntegratedFlux: True if users want to use time-integrated fluxes (integrated up to the investigated time point) in mole for analysis. Flase if users want to use instantaneous fluxes in mole/cm3/s at the investigated time point for analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "########## User Input ####################\n",
    "time_investigated = 0.000001 # hour\n",
    "reactant_name = 'ethane(1)'\n",
    "T = 1350 # K\n",
    "P = 100000 # Pa\n",
    "showIntegratedFlux = True \n",
    "##########################################\n",
    "\n",
    "\n",
    "timepoint_index = (np.abs(firstColDict['Time_(sec)']-time_investigated*3600)).argmin()\n",
    "timeList = firstColDict['Time_(sec)']\n",
    "Plist = firstColDict['Pressure_(bar)']\n",
    "Tlist = firstColDict['Temperature_(k)']\n",
    "Vlist = firstColDict['Volume_(cm3)']\n",
    "\n",
    "print \"Investigated time point is {0} secs\".format(timeList[timepoint_index])\n",
    "print \"At this moment, the reactant's mole fraction is {}.\".format(spc_mf_dict[reactant_name][timepoint_index])\n",
    "\n",
    "# Create moleFraction_dict at the investigated time point\n",
    "moleFraction_dict = {}\n",
    "for i, v in spc_mf_dict.items():\n",
    "    moleFraction_dict.update({i:v[timepoint_index]})\n",
    "\n",
    "# Create species_identifier_didct for the core species    \n",
    "species_identifier_dict = {}\n",
    "for i, s in enumerate(model.core.species):\n",
    "    species_identifier_dict.update({getSpeciesIdentifier(model.core.species[i]): s})\n",
    "\n",
    "# find the total consumption of the reactant in mol (positive value means consumption)\n",
    "reactantConsumption = 0 # \n",
    "for t in range(timepoint_index):\n",
    "    reactantROP0 = spc_total_rop_dict[reactant_name][1][t] # in mole/cm3*s\n",
    "    reactantROP1 = spc_total_rop_dict[reactant_name][1][t+1] # in mole/cm3*s\n",
    "    #integrate the flux using a trapezoidal rule\n",
    "    reactantConsumption += 0.5 * (reactantROP1 * Vlist[t+1] + reactantROP0 * Vlist[t]) * (timeList[t+1] - timeList[t])\n",
    "reactantConsumption = reactantConsumption * -1.0 # in mole\n",
    "\n",
    "if showIntegratedFlux:\n",
    "    print \"At this moment, the reactant's total molar consumption is {0:.3E} mole.\".format(reactantConsumption)\n",
    "\n",
    "# create DiGraph for this moment\n",
    "G = nx.DiGraph()\n",
    "for pair in graph_edges_dict:\n",
    "    node1 = getSpeciesIdentifier(pair[0])\n",
    "    node2 = getSpeciesIdentifier(pair[1])\n",
    "    e_rawdata = graph_edges_dict[pair]\n",
    "    total_flux = 0\n",
    "    for rxn in e_rawdata:\n",
    "        if showIntegratedFlux:\n",
    "            for t in range(timepoint_index):\n",
    "                rxnROP0 = e_rawdata[rxn][t] # in mole/cm3\n",
    "                rxnROP1 = e_rawdata[rxn][t+1] # in mole/cm3\n",
    "                #integrate the flux using a trapezoidal rule\n",
    "                total_flux += 0.5 * (rxnROP1 * Vlist[t+1] + rxnROP0 * Vlist[t]) * (timeList[t+1] - timeList[t]) # in mole\n",
    "        else:\n",
    "            total_flux += e_rawdata[rxn][timepoint_index] # in mole/cm3*s\n",
    "    if total_flux >= 0:\n",
    "        G.add_edge(node2, node1, {\"total_flux\":total_flux}) # in G, positive means production of node1\n",
    "    else: \n",
    "        G.add_edge(node1, node2, {\"total_flux\":-total_flux}) # in G, negative means consumption of node1   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [user input] Functionality 1: find most dominant species at the investigated time point\n",
    "\n",
    "Pleae specify\n",
    "\n",
    "- dominant_list: If you want to visualize the first top 10 dominant species, specify the list as [0, 9]. If you want to visualize from the 5th most dominant species to the 10th most dominant species, specify the list as [4, 9].\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "########## User Input ####################\n",
    "dominant_list = [0, 9] \n",
    "##########################################\n",
    "\n",
    "\n",
    "# sort out the moleFraction_dict in the descending order\n",
    "v = list(moleFraction_dict.values())\n",
    "k = list(moleFraction_dict.keys())\n",
    "index_list = [i[0] for i in sorted(enumerate(v), key=lambda x:x[1], reverse=True)]\n",
    "print '=============================================================='\n",
    "print 'Order \\t Species \\t Mole Frac. \\t MW (g/mol)'\n",
    "print '=============================================================='\n",
    "\n",
    "for order in range(dominant_list[0], dominant_list[-1]+1):\n",
    "    index_number = index_list[order]\n",
    "    species_dominant = k[index_number]\n",
    "    molecule_dominant = species_identifier_dict.get(species_dominant)\n",
    "    mole_fraction = v[index_number]\n",
    "    molecular_weight = molecule_dominant.molecule[0].getMolecularWeight()*1000\n",
    "    if len(species_dominant) < 6: # if the species label is too short, put more space to allign it\n",
    "        print '{0} \\t {1} \\t\\t {2} \\t {3:.0f}'.format(order, species_dominant, mole_fraction, molecular_weight)\n",
    "    else:\n",
    "        print '{0} \\t {1} \\t {2} \\t {3:.0f}'.format(order, species_dominant, mole_fraction, molecular_weight)\n",
    "    display(molecule_dominant)\n",
    "    print '--------------------------------------------------------------'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [user input] Functionality 2: find paths between two species\n",
    "\n",
    "Pleae specify\n",
    "\n",
    "- source species\n",
    "\n",
    "- target species"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "########## User Input ####################\n",
    "src_species = 'ethane(1)'\n",
    "tgt_species = 'C#C(25)'\n",
    "##########################################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## output1: sort out all the paths and most significant one"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "paths = list(nx.all_simple_paths(G, source=src_species, target=tgt_species, cutoff=5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "path_fluxes = []\n",
    "if showIntegratedFlux:\n",
    "    print '========================================================================================================'\n",
    "    print 'Path index \\t Integrated Flux (mole)    Path steps'\n",
    "else:\n",
    "    print '================================================'\n",
    "    print 'Path index \\t Flux (mole/cm3/s)    Path steps'\n",
    "print ''\n",
    "\n",
    "for i, path in enumerate(paths):\n",
    "    \n",
    "    path_steps = len(path) - 1\n",
    "    fluxes = [G[path[step]][path[step+1]]['total_flux'] for step in range(path_steps) ]\n",
    "    path_fluxes.append(min(fluxes))\n",
    "    if showIntegratedFlux:\n",
    "        print 'Path #{0}: \\t {1:.3E} \\t\\t  {2}'.format(i, min(fluxes), path)\n",
    "    else:\n",
    "        print 'Path #{0}: \\t {1:.3E} \\t {2}'.format(i, min(fluxes), path)\n",
    "sorted_path_fluxes = sorted(path_fluxes)\n",
    "\n",
    "print ''\n",
    "print '======================================================='\n",
    "print '\\t Most Significant Path '\n",
    "print '======================================================='\n",
    "if showIntegratedFlux:\n",
    "    print 'Path index \\t Integrated Flux (mole)'\n",
    "    print 'Path #{0}: \\t {1:.3E}'.format(path_fluxes.index(sorted_path_fluxes[-1]), sorted_path_fluxes[-1])\n",
    "else:\n",
    "    print 'Path index \\t Flux (mole/cm3/s)'\n",
    "    print 'Path #{0}: \\t {1:.3E} '.format(path_fluxes.index(sorted_path_fluxes[-1]), sorted_path_fluxes[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### [user input] output2: visualize the path of interest\n",
    "\n",
    "Please specify\n",
    "\n",
    "- the path index you want to visualize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "########## User Input ####################\n",
    "path_index_investigate = 0\n",
    "##########################################\n",
    "\n",
    "path = paths[path_index_investigate]\n",
    "path_steps = len(path) - 1\n",
    "print ''\n",
    "print '\\t Pathway Report '\n",
    "print '======================================'\n",
    "print 'The pathway you are intested in has {0} steps.'.format(len(path))\n",
    "for step in range(path_steps):\n",
    "    step_pair = (path[step], path[step+1])\n",
    "    h_abs_rxns = []\n",
    "    disp_rxns = []\n",
    "    Pdep_rxns = []\n",
    "    \n",
    "    print \"\"\n",
    "    if showIntegratedFlux:\n",
    "        print \"Step{0} \\t\\t\\t\\t Integrated Flux (mole)\".format(step+1)\n",
    "        print '***********************************************************************'\n",
    "        print \"{1} --> {2} \\t\\t {3:.3E}\".\\\n",
    "        format(step+1, step_pair[0], step_pair[1], G[step_pair[0]][step_pair[1]]['total_flux'])\n",
    "        print '***********************************************************************'\n",
    "    else:\n",
    "        print \"Step{0} \\t\\t\\t\\t Flux (mole/cm3/s)\".format(step+1)\n",
    "        print '*****************************************************'\n",
    "        print \"{1} --> {2} \\t\\t {3:.3E}\".\\\n",
    "        format(step+1, step_pair[0], step_pair[1], G[step_pair[0]][step_pair[1]]['total_flux'])\n",
    "        print '*****************************************************'\n",
    "    \n",
    "    if step_pair not in graph_edges_dict_simple:\n",
    "        step_pair = (step_pair[1], step_pair[0])\n",
    "    \n",
    "    print ''\n",
    "    for rxn in graph_edges_dict_simple[step_pair]:\n",
    "        if isinstance(rxn, PDepReaction):\n",
    "            Pdep_rxns.append(rxn)\n",
    "        elif rxn.family == \"H_Abstraction\":\n",
    "            h_abs_rxns.append(rxn)\n",
    "        elif rxn.family == \"Disproportionation\":\n",
    "            disp_rxns.append(rxn)\n",
    "        else:\n",
    "            print \"Example reaction: rxn#{0}\".format(rxn.index)\n",
    "            print ''\n",
    "            print str(rxn)\n",
    "            display(rxn)\n",
    "            print 'Forward rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, rxn.getRateCoefficient(T))\n",
    "            reverseRate = rxn.generateReverseRateCoefficient()\n",
    "            print 'Reverse rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, reverseRate.getRateCoefficient(T))\n",
    "            print ''\n",
    "    if len(h_abs_rxns) > 0: \n",
    "        \n",
    "        print \"Example H_Abstraction: rxn#{0} (1/{1} H_Abs): \".format(h_abs_rxns[0].index, len(h_abs_rxns)) \n",
    "        print ''\n",
    "        print str(h_abs_rxns[0])\n",
    "        display(h_abs_rxns[0])\n",
    "        print 'Forward rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, h_abs_rxns[0].getRateCoefficient(T))\n",
    "        reverseRate = h_abs_rxns[0].generateReverseRateCoefficient()\n",
    "        print 'Reverse rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, reverseRate.getRateCoefficient(T))\n",
    "        print ''\n",
    "    if len(disp_rxns) > 0: \n",
    "        \n",
    "        print \"Example Disproportionation: rxn#{0} (1/{1} Disp): \".format(disp_rxns[0].index, len(disp_rxns)) \n",
    "        print ''\n",
    "        print str(disp_rxns[0])\n",
    "        display(disp_rxns[0])\n",
    "        print 'Forward rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, disp_rxns[0].getRateCoefficient(T))\n",
    "        reverseRate = disp_rxns[0].generateReverseRateCoefficient()\n",
    "        print 'Reverse rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, reverseRate.getRateCoefficient(T))\n",
    "        print ''\n",
    "    if len(Pdep_rxns) > 0: \n",
    "        \n",
    "        print \"Example Pressure-dependent Rxn: rxn#{0} (1/{1} Pdep_rxns): \".format(Pdep_rxns[0].index, len(Pdep_rxns)) \n",
    "        print ''\n",
    "        print str(Pdep_rxns[0])\n",
    "        display(Pdep_rxns[0])\n",
    "        print 'Forward rate coefficient at {0} K and {1} Pa: {2:.2E} [cm3, mole, s]'.format(T, P, Pdep_rxns[0].getRateCoefficient(T, P))\n",
    "        reverseRate = Pdep_rxns[0].generateReverseRateCoefficient()\n",
    "        print 'Reverse rate coefficient at {0} K and {1} Pa: {2:.2E} [cm3, mole, s]'.format(T, P, reverseRate.getRateCoefficient(T, P))\n",
    "        print ''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [user input] Functionality 3: find paths from a particular species\n",
    "\n",
    "Please specify\n",
    "\n",
    "- source species\n",
    "\n",
    "- depth of pathway you want search, e.g., depth = 2 means search pathways 2 steps down from source species\n",
    "\n",
    "- path_top_list, e.g., [0, 3] means 1st step of the pathway takes the most significant branch and 2nd step would be the third most siginificant branch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "########## User Input ####################\n",
    "source = \"ethane(1)\"\n",
    "depth = 2\n",
    "path_top_list = [0, 0]\n",
    "##########################################\n",
    "\n",
    "\n",
    "current_node = source\n",
    "\n",
    "print ''\n",
    "print '\\t Pathway Report '\n",
    "print '======================================'\n",
    "print 'The pathway you are intested in has {0} steps.'.format(depth)\n",
    "for step in range(depth):\n",
    "    print \"\\n\"    \n",
    "    nextNode_flux_list = [(next_node, G[current_node][next_node]['total_flux']) for next_node in G[current_node]]\n",
    "    sorted_nextNode_flux_list = sorted(nextNode_flux_list, key=lambda tup: -tup[1])\n",
    "    \n",
    "    # choose the top one as next node\n",
    "    tup = sorted_nextNode_flux_list[path_top_list[step]]\n",
    "    next_node = tup[0]\n",
    "    step_flux = tup[1]\n",
    "    \n",
    "    print \"\"\n",
    "    if showIntegratedFlux:\n",
    "        print \"Step{0} \\t\\t\\t\\t Integrated Flux (mole)\".format(step+1)\n",
    "        print '************************************************************************'\n",
    "        print \"{0} --> {1} \\t\\t {2:.3E}\".\\\n",
    "        format(current_node, next_node, step_flux)\n",
    "        print '************************************************************************'\n",
    "    else:\n",
    "        print \"Step{0} \\t\\t\\t\\t Flux (mole/cm3/s)\".format(step+1)\n",
    "        print '****************************************************'\n",
    "        print \"{0} --> {1} \\t\\t {2:.3E}\".\\\n",
    "        format(current_node, next_node, step_flux)\n",
    "        print '****************************************************'\n",
    "\n",
    "    \n",
    "    step_pair = (current_node, next_node)\n",
    "    if step_pair not in graph_edges_dict_simple:\n",
    "        step_pair = (next_node, current_node)\n",
    "    \n",
    "    h_abs_rxns = []\n",
    "    disp_rxns = []\n",
    "    Pdep_rxns = []\n",
    "    for rxn in graph_edges_dict_simple[step_pair]:\n",
    "        if isinstance(rxn, PDepReaction):\n",
    "            Pdep_rxns.append(rxn)\n",
    "        elif rxn.family == \"H_Abstraction\":\n",
    "            h_abs_rxns.append(rxn)\n",
    "        elif rxn.family == \"Disproportionation\":\n",
    "            disp_rxns.append(rxn)\n",
    "        else:\n",
    "            print ''\n",
    "            print \"Example reaction: rxn#{0}\".format(rxn.index)\n",
    "            print ''\n",
    "            if showIntegratedFlux:\n",
    "            #integrate the flux using a trapezoidal rule\n",
    "                integrated_flux = 0\n",
    "                for t in range(timepoint_index):\n",
    "                    rxnROP0 = graph_edges_dict_simple[step_pair][rxn][t]\n",
    "                    rxnROP1 = graph_edges_dict_simple[step_pair][rxn][t+1]\n",
    "                    integrated_flux += 0.5 * (rxnROP1 * Vlist[t+1] + rxnROP0 * Vlist[t]) * (timeList[t+1] - timeList[t]) # in mole\n",
    "                print \"{0}: Integrated Flux = {1:.3E} mole\".format(str(rxn), integrated_flux)\n",
    "            else:\n",
    "                print \"{0}: Flux = {1:.3E} mole/cm3/s\".format(str(rxn), graph_edges_dict_simple[step_pair][rxn][timepoint_index])\n",
    "            display(rxn)\n",
    "            print 'Forward rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, rxn.getRateCoefficient(T))\n",
    "            reverseRate = rxn.generateReverseRateCoefficient()\n",
    "            print 'Reverse rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, reverseRate.getRateCoefficient(T))\n",
    "            print ''\n",
    "    if len(h_abs_rxns) > 0:\n",
    "        print ''\n",
    "        print \"Example H_Abstraction: rxn#{0}(1/{1} H_Abs)\".format(h_abs_rxns[0].index, len(h_abs_rxns))\n",
    "        print ''\n",
    "        print str(h_abs_rxns[0])\n",
    "        display(h_abs_rxns[0])\n",
    "        print 'Forward rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, h_abs_rxns[0].getRateCoefficient(T))\n",
    "        reverseRate = h_abs_rxns[0].generateReverseRateCoefficient()\n",
    "        print 'Reverse rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, reverseRate.getRateCoefficient(T))\n",
    "        print ''\n",
    "    if len(disp_rxns) > 0: \n",
    "        print ''\n",
    "        print \"Example Disproportionation: rxn#{0}(1/{1} Disp)\".format(disp_rxns[0].index, len(disp_rxns))\n",
    "        print ''\n",
    "        print str(disp_rxns[0])\n",
    "        display(disp_rxns[0])\n",
    "        print 'Forward rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, disp_rxns[0].getRateCoefficient(T))\n",
    "        reverseRate = disp_rxns[0].generateReverseRateCoefficient()\n",
    "        print 'Reverse rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, reverseRate.getRateCoefficient(T))\n",
    "        print ''\n",
    "    if len(Pdep_rxns) > 0: \n",
    "        print ''\n",
    "        print \"Example Pressure-dependent Rxn: rxn#{0} (1/{1} Pdep_rxns): \".format(Pdep_rxns[0].index, len(Pdep_rxns)) \n",
    "        print ''\n",
    "        print str(Pdep_rxns[0])\n",
    "        display(Pdep_rxns[0])\n",
    "        print 'Forward rate coefficient at {0} K and {1} Pa: {2:.2E} [cm3, mole, s]'.format(T, P, Pdep_rxns[0].getRateCoefficient(T, P))\n",
    "        reverseRate = Pdep_rxns[0].generateReverseRateCoefficient()\n",
    "        print 'Reverse rate coefficient at {0} K and {1} Pa: {2:.2E} [cm3, mole, s]'.format(T, P, reverseRate.getRateCoefficient(T, P))\n",
    "        print ''\n",
    "    \n",
    "    current_node = next_node"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### More detail for the H_Abstration and Disproportionation in last step of the pathway"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print ''\n",
    "print 'The step you see currently is between the two species {0}'.format(step_pair)\n",
    "print ''\n",
    "print 'Positive flux in below means positive production of 1st node of the pair.'\n",
    "print ''\n",
    "print 'This step have {0} H_Abstration and Disproportionation in total.'.format(len(h_abs_rxns + disp_rxns))\n",
    "\n",
    "flux_threshold = 1e-9\n",
    "\n",
    "print ''\n",
    "print 'TO avoid too much printout, the reactions shown below have fluxes above {0}.'.format(flux_threshold)\n",
    "\n",
    "total_flux = 0\n",
    "rxn_flux_tups = []\n",
    "for rxn in h_abs_rxns + disp_rxns:\n",
    "    if showIntegratedFlux:\n",
    "    #integrate the flux using a trapezoidal rule\n",
    "        flux = 0\n",
    "        for t in range(timepoint_index):\n",
    "            rxnROP0 = graph_edges_dict_simple[step_pair][rxn][t]\n",
    "            rxnROP1 = graph_edges_dict_simple[step_pair][rxn][t+1]\n",
    "            flux += 0.5 * (rxnROP1 * Vlist[t+1] + rxnROP0 * Vlist[t]) * (timeList[t+1] - timeList[t]) # in mole\n",
    "    else:\n",
    "        flux = graph_edges_dict_simple[step_pair][rxn][timepoint_index]\n",
    "    rxn_flux_tups.append((rxn, flux))\n",
    "    \n",
    "\n",
    "rxn_flux_tups = sorted(rxn_flux_tups, key=lambda tup: tup[1], reverse=False)\n",
    "for tup in rxn_flux_tups:\n",
    "    rxn = tup[0]\n",
    "    flux = tup[1]\n",
    "    if flux > flux_threshold:\n",
    "        total_flux += flux\n",
    "        print ''\n",
    "        print \"**********************************************************************************\"\n",
    "        if showIntegratedFlux:\n",
    "            print \"rxn#{0}: {1}: Integrated Flux = {2:.3E} mole \".format(rxn.index, str(rxn), flux) \n",
    "        else:\n",
    "            print \"rxn#{0}: {1}: Flux = {2:.3E} mole/cm3/s\".format(rxn.index, str(rxn), flux)\n",
    "        display(rxn) \n",
    "print \"***************************************\"\n",
    "print ''\n",
    "if showIntegratedFlux:\n",
    "    print \"TOTAL integrated flux from h_abs and disp is {0:.3E} mole.\".format(total_flux)\n",
    "else:\n",
    "    print \"TOTAL flux from h_abs and disp is {0:.3E} mole/cm3/s.\".format(total_flux)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [user input] Functionality 4: find paths towards a particular species\n",
    "\n",
    "Similar to functionality 2, please specify\n",
    "\n",
    "Please specify\n",
    "\n",
    "- target species\n",
    "\n",
    "- depth of pathway you want search, e.g., depth = 2 means search pathways 2 steps up beyond target species\n",
    "\n",
    "- path_top_list, e.g., [0, 3] means 1st step of the pathway takes the most significant branch and 2nd step would be the third most siginificant branch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "########## User Input ####################\n",
    "target = \"C#C(25)\"\n",
    "depth = 2\n",
    "path_top_list = [0, 0]\n",
    "##########################################\n",
    "\n",
    "\n",
    "current_node = target\n",
    "print ''\n",
    "print '\\t Pathway Report '\n",
    "print '======================================'\n",
    "print 'The pathway you are intested in has {0} steps.'.format(depth)\n",
    "print ''\n",
    "if showIntegratedFlux:\n",
    "    product_flux = 0\n",
    "    for t in range(timepoint_index):\n",
    "        productROP0 = spc_total_rop_dict[target][1][t]\n",
    "        productROP1 = spc_total_rop_dict[target][1][t+1]\n",
    "        product_flux += 0.5 * (productROP1 * Vlist[t+1] + productROP0 * Vlist[t]) * (timeList[t+1] - timeList[t]) # in mole\n",
    "    print \"total integrated product flux for {0} is {1:.3E} mole.\".format(target, product_flux)\n",
    "else:\n",
    "    print \"total product flux for {0} is {1:.3E} mole/cm3/s.\".format(target, spc_total_rop_dict[target][1][timepoint_index])\n",
    "for step in range(depth):\n",
    "    print \"\\n\"\n",
    "    prev_nodes = []\n",
    "    for node1 in G:\n",
    "        if current_node in G[node1]:\n",
    "            prev_nodes.append(node1)\n",
    "    prevNode_flux_list = [(prev_node, G[prev_node][current_node]['total_flux']) for prev_node in prev_nodes]\n",
    "    sorted_prevNode_flux_list = sorted(prevNode_flux_list, key=lambda tup: -tup[1])\n",
    "    \n",
    "    # choose the top one as next node\n",
    "    tup = sorted_prevNode_flux_list[path_top_list[step]]\n",
    "    prev_node = tup[0]\n",
    "    step_flux = tup[1]\n",
    "    \n",
    "    print \"\"\n",
    "    if showIntegratedFlux:\n",
    "        print \"Step{0} \\t\\t\\t\\t Integrated Flux (mole)\".format(step+1)\n",
    "        print '************************************************************************'\n",
    "        print \"{0} <-- {1} \\t\\t {2:.3E}\".\\\n",
    "        format(current_node, next_node, step_flux)\n",
    "        print '************************************************************************'\n",
    "    else:\n",
    "        print \"Step{0} \\t\\t\\t\\t Flux (mole/cm3/s)\".format(step+1)\n",
    "        print '****************************************************'\n",
    "        print \"{0} <-- {1} \\t\\t {2:.3E}\".\\\n",
    "        format(current_node, next_node, step_flux)\n",
    "        print '****************************************************'\n",
    "    \n",
    "    step_pair = (prev_node, current_node)\n",
    "    if step_pair not in graph_edges_dict_simple:\n",
    "        step_pair = (current_node, prev_node)\n",
    "    \n",
    "    h_abs_rxns = []\n",
    "    disp_rxns = []\n",
    "    Pdep_rxns = []\n",
    "    for rxn in graph_edges_dict_simple[step_pair]:\n",
    "        if isinstance(rxn, PDepReaction):\n",
    "            Pdep_rxns.append(rxn)\n",
    "        elif rxn.family == \"H_Abstraction\":\n",
    "            h_abs_rxns.append(rxn)\n",
    "        elif rxn.family == \"Disproportionation\":\n",
    "            disp_rxns.append(rxn)\n",
    "        else:\n",
    "            print ''\n",
    "            print \"Example reaction: rxn#{0}\".format(rxn.index)\n",
    "            print ''\n",
    "            if showIntegratedFlux:\n",
    "            #integrate the flux using a trapezoidal rule\n",
    "                integrated_flux = 0\n",
    "                for t in range(timepoint_index):\n",
    "                    rxnROP0 = graph_edges_dict_simple[step_pair][rxn][t]\n",
    "                    rxnROP1 = graph_edges_dict_simple[step_pair][rxn][t+1]\n",
    "                    integrated_flux += 0.5 * (rxnROP1 * Vlist[t+1] + rxnROP0 * Vlist[t]) * (timeList[t+1] - timeList[t]) # in mole\n",
    "                print \"{0}: Integrated Flux = {1:.3E} mole\".format(str(rxn), integrated_flux)\n",
    "            else:\n",
    "                print \"{0}: Flux = {1:.3E} mole/cm3/s\".format(str(rxn), graph_edges_dict_simple[step_pair][rxn][timepoint_index])\n",
    "            display(rxn)\n",
    "            print 'Forward rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, rxn.getRateCoefficient(T))\n",
    "            reverseRate = rxn.generateReverseRateCoefficient()\n",
    "            print 'Reverse rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, reverseRate.getRateCoefficient(T))\n",
    "            \n",
    "    if len(h_abs_rxns) > 0:\n",
    "        print ''\n",
    "        print \"Example H_Abstraction: rxn#{0}(1/{1} H_Abs)\".format(h_abs_rxns[0].index, len(h_abs_rxns))\n",
    "        print ''\n",
    "        print str(h_abs_rxns[0])\n",
    "        display(h_abs_rxns[0])\n",
    "        print 'Forward rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, h_abs_rxns[0].getRateCoefficient(T))\n",
    "        reverseRate = h_abs_rxns[0].generateReverseRateCoefficient()\n",
    "        print 'Reverse rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, reverseRate.getRateCoefficient(T))\n",
    "        print ''\n",
    "    if len(disp_rxns) > 0: \n",
    "        print ''\n",
    "        print \"Example Disproportionation: rxn#{0}(1/{1} Disp)\".format(disp_rxns[0].index, len(disp_rxns))\n",
    "        print ''\n",
    "        print str(disp_rxns[0])\n",
    "        display(disp_rxns[0])\n",
    "        print 'Forward rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, disp_rxns[0].getRateCoefficient(T))\n",
    "        reverseRate = disp_rxns[0].generateReverseRateCoefficient()\n",
    "        print 'Reverse rate coefficient at {0} K: {1:.2E} [cm3, mole, s]'.format(T, reverseRate.getRateCoefficient(T))\n",
    "        print ''\n",
    "    if len(Pdep_rxns) > 0: \n",
    "        print ''\n",
    "        print \"Example Pressure-dependent Rxn: rxn#{0} (1/{1} Pdep_rxns): \".format(Pdep_rxns[0].index, len(Pdep_rxns)) \n",
    "        print ''\n",
    "        print str(Pdep_rxns[0])\n",
    "        display(Pdep_rxns[0])\n",
    "        print 'Forward rate coefficient at {0} K and {1} Pa: {2:.2E} [cm3, mole, s]'.format(T, P, Pdep_rxns[0].getRateCoefficient(T, P))\n",
    "        reverseRate = Pdep_rxns[0].generateReverseRateCoefficient()\n",
    "        print 'Reverse rate coefficient at {0} K and {1} Pa: {2:.2E} [cm3, mole, s]'.format(T, P, reverseRate.getRateCoefficient(T, P))\n",
    "        print ''\n",
    "    \n",
    "    current_node = prev_node"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### More detail for the H_Abstration and Disproportionation in last step of the pathway"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print ''\n",
    "print 'The step you see currently is between the two species {0}'.format(step_pair)\n",
    "print ''\n",
    "print 'Positive flux in below means positive production of 1st node of the pair.'\n",
    "print ''\n",
    "print 'This step have {0} H_Abstration and Disproportionation in total.'.format(len(h_abs_rxns + disp_rxns))\n",
    "\n",
    "flux_threshold = 1e-9\n",
    "\n",
    "print ''\n",
    "print 'TO avoid too much printout, the reactions shown below have fluxes above {0}.'.format(flux_threshold)\n",
    "\n",
    "total_flux = 0\n",
    "rxn_flux_tups = []\n",
    "for rxn in h_abs_rxns + disp_rxns:\n",
    "    if showIntegratedFlux:\n",
    "    #integrate the flux using a trapezoidal rule\n",
    "        flux = 0\n",
    "        for t in range(timepoint_index):\n",
    "            rxnROP0 = graph_edges_dict_simple[step_pair][rxn][t]\n",
    "            rxnROP1 = graph_edges_dict_simple[step_pair][rxn][t+1]\n",
    "            flux += 0.5 * (rxnROP1 * Vlist[t+1] + rxnROP0 * Vlist[t]) * (timeList[t+1] - timeList[t]) # in mole\n",
    "    else:\n",
    "        flux = graph_edges_dict_simple[step_pair][rxn][timepoint_index]\n",
    "    rxn_flux_tups.append((rxn, flux))\n",
    "\n",
    "rxn_flux_tups = sorted(rxn_flux_tups, key=lambda tup: tup[1], reverse=False)\n",
    "for tup in rxn_flux_tups:\n",
    "    rxn = tup[0]\n",
    "    flux = tup[1]\n",
    "    if flux > flux_threshold:\n",
    "        total_flux += flux\n",
    "        print ''\n",
    "        print \"**********************************************************************************\"\n",
    "        if showIntegratedFlux:\n",
    "            print \"rxn#{0}: {1}: Integrated Flux = {2:.3E} mole \".format(rxn.index, str(rxn), flux) \n",
    "        else:\n",
    "            print \"rxn#{0}: {1}: Flux = {2:.3E} mole/cm3/s\".format(rxn.index, str(rxn), flux)\n",
    "        display(rxn) \n",
    "print \"***************************************\"\n",
    "print ''\n",
    "if showIntegratedFlux:\n",
    "    print \"TOTAL integrated flux from h_abs and disp is {0:.3E} mole.\".format(total_flux)\n",
    "else:\n",
    "    print \"TOTAL flux from h_abs and disp is {0:.3E} mole/cm3/s.\".format(total_flux)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
